iv3 = bonds_GDP/deposits_GDP
iv4 = NA
mm = cbind( mm0=residual, mm_VIX=residual*iv1, mm_deposits_spread=residual*iv2, mm_bonds_to_deposits=residual*iv3 )
}
if(shut_off_VIX){
column_list = 1:ncol(mm)
mm = mm[ , column_list[-2] ]  # Remove the second moment with VIX
}
if(!prediction){
return(mm*weights)
}else{
liq_prem_predicted_Only_FFR =   deposits_spread * mean( (lambda)*( bonds_GDP/deposits_GDP )^(rho-1) )
liq_prem_predicted_Only_Supply = mean(deposits_spread) * (mean(lambda)*( bonds_GDP/mean(deposits_GDP) )^(rho-1) )
liq_prem_predicted_Only_VIX = mean(deposits_spread) * (lambda* mean(( bonds_GDP/deposits_GDP )^(rho-1)) )
liq_prem_predicted_Only_Supply_and_VIX = mean(deposits_spread) * (lambda* ( bonds_GDP/mean(deposits_GDP) )^(rho-1) )
liq_prem_predicted_Only_Supply_and_FFR = deposits_spread * (mean(lambda)* ( bonds_GDP/mean(deposits_GDP) )^(rho-1) )
liq_prem_predicted_Only_VIX_and_FFR = deposits_spread * ((lambda)* mean(( bonds_GDP/deposits_GDP )^(rho-1)) )
LiqPrem_residual= LiqPrem-liq_prem_predicted
deposits_spread_predicted = LiqPrem / (lambda*( bonds_GDP/deposits_GDP )^(rho-1) )
return( cbind(data.table( residual, LiqPrem, liq_prem_predicted, deposits_spread, deposits_spread_predicted, LiqPrem_residual, iv1, iv2, iv3, iv4,  liq_prem_predicted_Only_FFR,liq_prem_predicted_Only_Supply,liq_prem_predicted_Only_VIX,liq_prem_predicted_Only_Supply_and_VIX,liq_prem_predicted_Only_Supply_and_FFR,liq_prem_predicted_Only_VIX_and_FFR), as.data.frame(mm)  ) )
}
}
# ===================== TABLE 4: Core GMM Estimations =====================
beta_init = c( 0.3, 0.01 )
DT = TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX, DebtToGDP_raw ) ]
reg = gmm_customized(g_single, DT, beta_init, optfct = "nlminb", lower=c(0, 0), upper=c(30,1) )
result = g_single(reg$coefficients, DT, prediction = TRUE)
R2_vec = R2_fun( result$liq_prem_predicted, result$LiqPrem )
core_GMM_fun <- function(TBs, use_debt_GDP_IV=FALSE, shut_off_VIX=FALSE, beta_init=c(0.4,0.01), lags=4, use_AR1=FALSE, use_alternative_residual=FALSE, type="iterative"){
N = length(TBs)
R2_vec = rep(0, N)
p_value_vec = rep(0, N)
regs = list()
for(iter in 1:N){
DT = remove_missing(as.data.frame(TBs[[iter]]))
reg = gmm_customized( function(beta,X) {g_single(beta, X, use_debt_GDP_IV = use_debt_GDP_IV, shut_off_VIX=shut_off_VIX, use_AR1 = use_AR1, use_alternative_residual = use_alternative_residual)} , DT, beta_init, optfct = "nlminb", bw=lags, type=type )
result = g_single(reg$coefficients, DT, prediction = TRUE, use_debt_GDP_IV = use_debt_GDP_IV, shut_off_VIX=shut_off_VIX, use_AR1 = use_AR1, use_alternative_residual = use_alternative_residual)
R2_vec[iter] = R2_fun( result$liq_prem_predicted, result$LiqPrem )
p_value_vec[iter] = round( as.numeric(summary(reg)$stest$test[2])  , digits=3  )
regs[[iter]] = reg
}
return( list( regs=regs, R2_vec=R2_vec, p_value_vec=p_value_vec )  )
}
# Regression Settings
TBs = list( TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ],
TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ]
)
result = core_GMM_fun(TBs)
result_baseline = result
result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE)
result_with_IV_baseline = result_with_IV
result_core_list = c(result$regs, result_with_IV$regs)
reg_baseline = result_with_IV$regs[[1]]
#  ***** Main GMM Table *****
# Latex table
stargazer( c(result$regs, result_with_IV$regs), omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt",
covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines =
list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ),
c( "Total Treasury IV?", "No", "No", "Yes", "Yes" ),
c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
column.labels  = rep(c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),2) ,
dep.var.caption = "Measurement of B/D" )
# Print output
stargazer( c(result$regs, result_with_IV$regs) , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt",
covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines =
list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ),
c( "IV of $ \\frac{\\text{total debt}}{\\text{GDP}} $?", "No", "No", "Yes", "Yes" ),
c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
column.labels  = rep(c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),2) ,
dep.var.caption = "Measurement of B/D", type="text" )
# ==== RFS Revision Table: AR1 Coefficients ====
result = core_GMM_fun(TBs, use_AR1 = TRUE, beta_init=c(0.4,0.01,0.6))
result_with_IV = core_GMM_fun(TBs,use_AR1 = TRUE, use_debt_GDP_IV = TRUE, beta_init=c(0.4,0.01,0.6))
stargazer( c(result$regs, result_with_IV$regs) , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt",
covariate.labels = c( "$\\rho$", "$\\beta$", "\\zeta" ), add.lines =
list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ),
c( "IV of $ \\frac{\\text{total debt}}{\\text{GDP}} $?", "No", "No", "Yes", "Yes" ),
c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
column.labels  = rep(c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),2) ,
dep.var.caption = "Measurement of B/D", type="text" )
# Autocorrelations
result_ACF = g_single( c(result$regs[[1]]$coefficients[1:2],0) , TBs[[1]], prediction = TRUE, use_AR1 = FALSE)
plot(result_ACF$residual)
ggAcf(result_ACF$residual) + ggtitle("")
ggsave(width = 6, height = 4, dpi = 300, filename = paste0(FigurePath,"ACF_plot_of_residual.pdf"))
durbinWatsonTest(result_ACF$residual, max.lag=2)
# # GMM with multiplicative residual (divided by quantity ratio and VIX)
# # ==== RFS Revision Table: GMM multiplicative residuals ====
# result = core_GMM_fun(TBs, use_alternative_residual = TRUE)
# result_with_IV = core_GMM_fun(TBs, use_alternative_residual = TRUE, use_debt_GDP_IV = TRUE)
# GMM without VIX
# ===== RFS Revision Table: GMM without VIX =====
result = core_GMM_fun(TBs, lags=4, shut_off_VIX=TRUE)
result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE, lags=4, shut_off_VIX=TRUE)
# ===================== TABLE 5: Explanatory Power of Different Components =====================
results = g_single( reg_baseline$coefficients, TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ], prediction=TRUE )
R2s = abs(c( R2_fun( results$liq_prem_predicted_Only_Supply, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_FFR, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_VIX, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_Supply_and_FFR, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_VIX_and_FFR, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted, DT$LiqPrem  )  ))
print( paste0("R2 = ", paste0(round(100*R2s,digits = 1),collapse="%, "), "%"  ), digits = 3 )
# ===================== Figure: Basic GMM Illustrations =====================
result = g_single( reg$coefficients, DT[, c(1:5), with=FALSE ], prediction = TRUE )
# Panel A: Prediction on Liquidity premium
pdf( paste0(FigurePath,"Model_and_Data.pdf"), width = 5, height = 4 )
par(mar = c(2, 3.5, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, result[ , list(LiqPrem, liq_prem_predicted) ]  , ylab="Liquidity Premium (%)", LegendNames = c("Data", "Model Predictions") , LegendPosition = "topright" )
abline( v=as.Date(paste0(WWII_periods[c(1,length(WWII_periods))],"-01-01")), lty=2 )
text( as.numeric( as.Date( paste0(round(mean(WWII_periods)), "-01-01") ) ) , 2, labels = "WWII"   )
dev.off()
# Panel B: Prediction on FFR
pdf( paste0(FigurePath,"Model_and_Data2.pdf"), width = 5, height = 4 )
par(mar = c(2, 3.5, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, result[ , list(deposits_spread/deposits_spread_sensitivity_to_FFR, deposits_spread_predicted/deposits_spread_sensitivity_to_FFR) ]  , ylab="Federal Funds Rate (%)", LegendNames = c("Data", "Model Predictions") , LegendPosition = "topright" )
abline( v=as.Date(paste0(WWII_periods[c(1,length(WWII_periods))],"-01-01")), lty=2 )
text( as.numeric( as.Date( paste0(round(mean(WWII_periods)), "-01-01") ) ) , 20, labels = "WWII"   )
dev.off()
# ===================== Figure: Liquidity Share of Treasuries (Plot the beta*VIX) ==============
# Note: check the ratio
Ratios = reg_baseline$coefficients[2] * TB_preGMM$VIX
lambdas_measured = Ratios/(1+Ratios)
mean(lambdas_measured)
library(latex2exp)
pdf( paste0(FigurePath,"Liquidity_Share_of_Treasuries.pdf"), width = 7, height = 4.5 )
par(mar = c(2, 3.5, 0.5, 0.5), mfrow=c(1,1))
lambdas_smooothed = rollmean(lambdas_measured, 120, fill = "extend")
MyPlot( TB_preGMM$YearMon, data.frame(lambdas_measured,lambdas_smooothed ), ylim=c(0,0.5), ylab = "Liquidity Share of Treasuries",  LegendNames = c(expression(paste("Measured ", lambda[t])),"Moving Average (10 years)") )
dev.off()
mean(lambdas_smooothed)
# ===================== TABLE 6: Subsample Analyses =====================
cutoff_date = min( TB_preGMM$YearMon[!is.na(TB_preGMM$Deposit_Spread_raw)] )
TBs = list( TB_preGMM[ !is.element(Year, WWII_periods)  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ],
TB_preGMM[ !is.element(Year, WWII_periods)  , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ],
TB_preGMM[ Year>max(WWII_periods) , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ],
TB_preGMM[ Year>max(WWII_periods) , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ],
TB_preGMM[ Year>1980 , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ],
TB_preGMM[ Year>1980 , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ]
)
result = core_GMM_fun(TBs)
result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE)
stargazer( c(result_with_IV$regs) , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt",
covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines =
list( c("p-value of J-test",  c(result_with_IV$p_value_vec) ),
c( "Variations explained", paste0(round( c(result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
column.labels  = rep( c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),3) ,
dep.var.caption = "Measurement of B/D" )
# ================ **SECTION 3** GMM Algorithm for the Nested Model =============================
g_inner <- function( beta, X_matrix, prediction=FALSE, use_debt_GDP_IV=FALSE ){
X_matrix = as.matrix(X_matrix)
eta = beta[1]
beta_mu = beta[2]
# The structure of X_matrix = [  liquidity_premium (P2_CP_Tbill), Tsy/GDP, VIX, KVJ_net_deposits(non_bank_debt_GDP),  LiqPrem_nonbank=P2_P1_spread ,   DebtToGDP_raw,  credit_spread, FFR ]
LiqPrem_Tsy = X_matrix[,1]
Tsy_GDP = X_matrix[,2]
VIX = X_matrix[,3]
non_bank_debt_GDP = X_matrix[,4]
LiqPrem_non_bank = X_matrix[,5]
total_debt_to_GDP = X_matrix[,6]
FFR = X_matrix[,7]
LiqPrem_Tsy_predicted = beta_mu * VIX*(Tsy_GDP/non_bank_debt_GDP)^(eta-1)*LiqPrem_non_bank
residual_inner = LiqPrem_Tsy - LiqPrem_Tsy_predicted
iv1 = FFR
iv2 = VIX
if(use_debt_GDP_IV){
iv3 = (total_debt_to_GDP)
iv4 = (total_debt_to_GDP)^2
mm = cbind( residual_inner, residual_inner*iv1, residual_inner*iv2, residual_inner*iv3, residual_inner*iv4)
}else{
iv3 = Tsy_GDP/non_bank_debt_GDP
mm = cbind( residual_inner, residual_inner*iv1, residual_inner*iv2, residual_inner*iv3)
}
if(!prediction){
return(mm)
}else{
return( data.table( LiqPrem_Tsy, LiqPrem_Tsy_predicted) )
}
}
GMMs_inner_fun <- function(TBs, use_debt_GDP_IV=FALSE, rho_limit=2, remove_negative_liq_prem = FALSE, beta_init=c(0.5,0.1), type="twoStep"){
N = length(TBs)
R2_vec = rep(0, N)
p_value_vec = rep(0, N)
regs = list()
for(iter in 1:N){
DT = remove_missing(as.data.frame(TBs[[iter]]))
if(remove_negative_liq_prem){
negative_index = DT$LiqPrem<0 | DT$LiqPrem_nonbank<0
print( paste0( sum(negative_index), " observations removed due to negative liquidity premium values!" ) )
DT = DT[!negative_index,]
}
reg = gmm_customized( function(beta,X) {g_inner(beta, X, use_debt_GDP_IV = use_debt_GDP_IV)} , DT, beta_init, optfct = "nlminb", lower=c(0, 0), upper=c(rho_limit,1), type = type )
result = g_inner(reg$coefficients, DT, prediction = TRUE, use_debt_GDP_IV = use_debt_GDP_IV)
R2_vec[iter] = cor( result$LiqPrem_Tsy_predicted, result$LiqPrem_Tsy )^2
p_value_vec[iter] = round( summary(reg)$stest$test[2], digits=3  )
regs[[iter]] = reg
}
return( list( regs=regs, R2_vec=R2_vec, p_value_vec=p_value_vec )  )
}
sum(  TB_preGMM$P2CP3M -TB_preGMM$MMF_rate<0, na.rm=TRUE )
par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( AACP3M, P2CP3M ) ], LegendNames = c( "AA CP", "P2 CP" ), NA_elimination = FALSE )
pdf( paste0(FigurePath, "comparison_liq_prem.pdf"), width=8, height =4 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,2))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list(FFR - MMF_rate, Repo3M - MMF_rate, AACP3M - MMF_rate, P2CP3M - MMF_rate  ) ], LegendNames = c( "FFR-MMF rate", "Repo 3M-MMF rate", "AA CP-MMF rate", "P2 CP-MMF rate" ) )
abline(h=0)
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list(FFR - AACP3M, Repo3M - AACP3M, P2CP3M - AACP3M ) ], LegendNames = c( "FFR-AA CP rate", "Repo 3M-AA CP rate", "P2 CP-AA CP rate" ) )
abline(h=0)
dev.off()
pdf( paste0(FigurePath, "comparison_non_bank_volume.pdf"), width=5, height =4 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( MMF_to_GDP, KVJ_fin_sector_debt_to_GDP - Deposits_GDP ) ], LegendNames = c( "MMF/GDP", "(KVJ fin debt - deposits)/GDP" )  )
abline(h=0)
dev.off()
# Regression Settings
TB_reg = TB_preGMM[Year>=1976]
# Fed data: include P1 CP
par(mfrow=c(1,1))
MyPlot( TB_reg$YearMon, TB_reg[ , list( P2_CP_Tbill, P2_P1_spread,  P2_CP_Repo) ], LegendNames = c( "P2 CP - Tsy", "P2 CP - AA CP", "P2 CP - Repo" ) )
# Plot the relationship
par(mar = c(4, 4, 3, 3), mfrow=c(1,2))
plot( log(TB_reg$DebtToGDP/TB_reg$KVJ_net_deposits), log(TB_reg$P2_CP_Tbill/TB_reg$P2_P1_spread), xlab="log(Net Tsy/Net Non-bank Deposits)", ylab="log(P2CP-Tbill/P2-P1)" , mgp = c(2.2, 0.8, 0))
plot( log(TB_reg$DebtToGDP/TB_reg$KVJ_net_deposits), log(TB_reg$P2_CP_Tbill/TB_reg$P2_CP_Repo), xlab="log(Net Tsy/Net Non-bank Deposits)", ylab="log(P2CP-Tbill/P2-Repo)" , mgp = c(2.2, 0.8, 0) )
# Plot the ratios
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( P2_CP_Tbill/P2_P1_spread ) ], LegendNames = c( "P2 CP-Tbill/P2-P1 spread" ) )
pdf( paste0(FigurePath, "/CP_measures.pdf"), width=5, height = 4)
par(mar = c(3.5, 3.5, 1.0, 1.0), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( P2_CP_Tbill,  P2_P1_spread, P2_CP_Repo ) ], LegendNames = c( "P2 CP-Tsy 3M spread", "P2 CP-P1 CP 3M spread", "P2 CP-Repo 3M spread" ) )
abline(h=0)
dev.off()
# Do a log regression
regs = list()
regs[[1]] = lm(  log(P2_CP_Tbill/P2_P1_spread) ~ log(DebtToGDP/KVJ_net_deposits) + log(VIX), data=TB_reg  )
regs[[2]] = lm(  log(P2_CP_Tbill/P2_CP_Repo) ~ log(DebtToGDP/KVJ_net_deposits) + log(VIX), data=TB_reg  )
regs[[3]] = ivreg( log(P2_CP_Tbill/P2_P1_spread) ~ log(DebtToGDP/KVJ_net_deposits) + log(VIX) | log(P2_P1_spread) + DebtToGDP_raw + DebtToGDP_raw^2 + log(VIX) , data=TB_reg )
regs[[4]] = ivreg( log(P2_CP_Tbill/P2_CP_Repo) ~ log(DebtToGDP/KVJ_net_deposits) + log(VIX) | log(P2_P1_spread) + DebtToGDP_raw + DebtToGDP_raw^2 + log(VIX) , data=TB_reg )
stargazer(regs, se=RobustSE(regs),  star.cutoffs = c(0,0,0), type="text", omit.stat = c("adj.rsq","F", "ser"), dep.var.labels.include = FALSE, column.labels = rep(c("P2-P1", "P2-Repo"),2), dep.var.caption = "Dep var: log(Tsy liq prem/non-bank liq prem)" )
regs = list()
regs[[1]] = lm(  log(P2_CP_Tbill/P2_P1_spread) ~ log(DebtToGDP/KVJ_net_deposits) , data=TB_reg  )
regs[[2]] = lm(  log(P2_CP_Tbill/P2_CP_Repo) ~ log(DebtToGDP/KVJ_net_deposits) , data=TB_reg  )
regs[[3]] = ivreg( log(P2_CP_Tbill/P2_P1_spread) ~ log(DebtToGDP/KVJ_net_deposits)  | log(P2_P1_spread) + DebtToGDP_raw + DebtToGDP_raw^2 + log(VIX) , data=TB_reg )
regs[[4]] = ivreg( log(P2_CP_Tbill/P2_CP_Repo) ~ log(DebtToGDP/KVJ_net_deposits) | log(P2_P1_spread) + DebtToGDP_raw + DebtToGDP_raw^2 + log(VIX) , data=TB_reg )
stargazer(regs, se=RobustSE(regs),  star.cutoffs = c(0,0,0), type="text", omit.stat = c("adj.rsq","F", "ser"), dep.var.labels.include = FALSE, column.labels = rep(c("P2-P1", "P2-Repo"),2), dep.var.caption = "Dep var: log(Tsy liq prem/non-bank liq prem)" )
# ===================== Table 7: GMM Estimation of Inner Layer of the Nested System =============================
TB_reg = TB_preGMM[Year>=1970]
TBs = list( TB_reg[  , list( P2_CP_Tbill,    DebtToGDP , VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_P1_spread ,   DebtToGDP_raw,  credit_spread, FFR) ],
TB_reg[  , list( P2_CP_Tbill,    DebtToGDP, VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_CP_MMF,     DebtToGDP_raw,  credit_spread, FFR) ]
)
GMM_fun = GMMs_inner_fun
rho_upper = 1
result = GMM_fun(TBs, rho_limit = rho_upper, beta_init = c(0.5,0.1), type="twoStep")
result_with_IV = GMM_fun(TBs,use_debt_GDP_IV = TRUE, rho_limit = rho_upper,beta_init = c(0.5,0.1), type="iterative")
column_labels = rep(c("P2-P1 CP", "P2CP-MMF"),2)
lag_input = 4
stargazer( c(result$regs, result_with_IV$regs) ,  se=RobustSE(c(result$regs, result_with_IV$regs),lag_input = lag_input),
add.lines = list( c("Variation Explained",  paste0(round(100*c( result$R2_vec, result_with_IV$R2_vec ),digits=1),"\\%") ),
c( "Total Treasury IV?", rep("No",length(result$regs)), rep("Yes",length(result_with_IV$regs)) ) ),
dep.var.labels = paste0("GMM (NeweyWest with lag=", lag_input,")"),
dep.var.caption=c(""), star.cutoffs = c(0,0,0), covariate.labels= c("$\\eta$", "$\\beta_\\mu$"), column.labels = column_labels  )
# Use P2P1 spread as the non-bank liquidity premium measure
result_inner_list = c( result$regs, result_with_IV$regs )
result_inner = result_with_IV$regs[[1]]
eta = result_inner$coefficients[1]
beta_mu = result_inner$coefficients[2]
TB_preGMM[ , mu :=beta_mu*VIX/(1+beta_mu*VIX)  ]
TB_preGMM[ , bonds_GDP_inner :=  ( (1-mu)*KVJ_net_deposits^eta + mu*DebtToGDP^eta )^(1/eta)  ]
TB_preGMM[ , LiqPrem_bond_inner := ( (1-mu)^(-1/(eta-1))*P2_P1_spread^(eta/(eta-1)) + mu^(-1/(eta-1))*P2_CP_Tbill^(eta/(eta-1)) )^((eta-1)/eta)  ]
par(mar = c(3.5, 3.5, 1.0, 1.0), mfrow=c(1,2))
MyPlot(TB_preGMM$YearMon,  TB_preGMM[ , list( bonds_GDP_inner,DebtToGDP ) ], LegendNames = c("Constructed Inner Bond/GDP", "Net Tsy/GDP")  )
MyPlot(TB_preGMM$YearMon,  TB_preGMM[ , list( LiqPrem_bond_inner,P2_CP_Tbill ) ], LegendNames = c("Constructed Inner Bond LiqPrem", "Tsy Liq Prem")  )
MyPlot( TB$YearMon,  TB[ , list(Deposits_GDP +  DebtToGDP + KVJ_net_deposits  ) ] )
MyPlot( TB[ Year>1950]$YearMon,  TB[ Year>1950 , list(standardize(FFR), standardize(Deposits_GDP)  ) ], LegendNames = c("FFR", "Deposits/GDP"), ylab="Standardized values"  )
MyPlot( TB[ Year>1950]$YearMon,  TB[ Year>1950 , list(standardize(FFR), standardize(DebtToGDP )  ) ], LegendNames = c("FFR", "Net Tsy/GDP"), ylab="Standardized values" )
MyPlot( TB[ Year>1950]$YearMon,  TB[ Year>1950 , list(standardize(FFR), standardize(KVJ_net_deposits )  ) ], LegendNames = c("FFR", "Net Tsy/GDP"), ylab="Standardized values" )
# ===================== Table 8: GMM Estimation of the Whole Nested System =============================
g_nested <- function(beta,X_matrix, weights=1, prediction=FALSE, total_Tsy_IV=TRUE, use_VIX_moment=TRUE, alternative_residual=FALSE)
{
# The columns of X_matrix are: (1) Tsy liquidity premium, (2) deposits spread, (3) bonds/GDP, (4) deposits/GDP, (5) VIX, (6) MMF/GDP,  (7) LiqPrem_MMF,  (8) Total Debt/GDP, (9) FFR
X_matrix = as.data.frame(X_matrix)
rho = beta[1]
beta_lambda = beta[2]
eta = beta[3]
beta_mu = beta[4]
# The structure of X_matrix = [  liquidity_premium, Deposit_spread, Tsy/GDP, deposit/GDP, VIX,  nonbank_deposits_GDP,  LiqPrem_nonbank, total_debt_to_GDP,  FFR]
LiqPrem_Tsy = X_matrix[,1]
deposits_spread = X_matrix[,2]
Tsy_GDP = X_matrix[,3]
deposits_GDP = X_matrix[,4]
VIX = X_matrix[,5]
nonbank_deposits_GDP = X_matrix[,6]
LiqPrem_nonbank = X_matrix[,7]
total_debt_to_GDP = X_matrix[,8]
FFR = X_matrix[,9]
mu_t = beta_mu*VIX/(1+beta_mu*VIX)  # Important: remmeber how mu_t is approxed
LiqPrem_Tsy_predicted = mu_t/(1-mu_t)*(Tsy_GDP/nonbank_deposits_GDP)^(eta-1)*LiqPrem_nonbank
residual_inner = LiqPrem_Tsy - LiqPrem_Tsy_predicted
if(alternative_residual){
residual_inner = LiqPrem_Tsy/(mu_t/(1-mu_t)*(Tsy_GDP/nonbank_deposits_GDP)^(eta-1)) - LiqPrem_nonbank
}
LiqPrem_bonds = ( (1-mu_t)^(-1/(eta-1))*LiqPrem_nonbank^(eta/(eta-1)) + mu_t^(-1/(eta-1))*LiqPrem_Tsy^(eta/(eta-1)) )^((eta-1)/eta)
bonds_GDP = ( (1-mu_t)*nonbank_deposits_GDP^eta + mu_t*Tsy_GDP^eta )^(1/eta)
LiqPrem_bonds_predicted = deposits_spread*beta_lambda*VIX*( bonds_GDP/deposits_GDP )^(rho-1)
residual_outter = LiqPrem_bonds - LiqPrem_bonds_predicted
if(alternative_residual){
residual_outter = LiqPrem_bonds/(beta_lambda*VIX*( bonds_GDP/deposits_GDP )^(rho-1)) -  deposits_spread
}
iv_FFR = FFR
iv4 = total_debt_to_GDP
iv5 = (total_debt_to_GDP)^2
iv_q_inner = Tsy_GDP/nonbank_deposits_GDP
iv_q_outer = bonds_GDP/deposits_GDP
if(total_Tsy_IV){
if(use_VIX_moment){
mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank, residual_inner*VIX, residual_inner*iv_q_inner,
residual_outter, residual_outter*iv_FFR, residual_outter*VIX, residual_outter*iv4, residual_outter*iv5 )
}else{
mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank,  residual_inner*iv_q_inner,
residual_outter, residual_outter*iv_FFR, residual_outter*iv4, residual_outter*iv5 )
}
}else{
if(use_VIX_moment){
mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank, residual_inner*VIX, residual_inner*iv_q_inner,
residual_outter, residual_outter*iv_FFR, residual_outter*VIX, residual_outter*iv_q_outer )
}else{
mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank, residual_inner*iv_q_inner,
residual_outter, residual_outter*iv_FFR,  residual_outter*iv_q_outer )
}
}
if(!prediction){
return(mm*weights)
}else{
return( data.table( LiqPrem_bonds, LiqPrem_bonds_predicted, LiqPrem_Tsy, LiqPrem_Tsy_predicted) )
}
}
coef_estimations = regs[[1]]$coefficients
R2_fun_nested <- function( coef_estimations, DT, total_Tsy_IV=TRUE ){
result = g_nested(coef_estimations, DT, prediction=TRUE, total_Tsy_IV=total_Tsy_IV)
R2_Tsy = R2_fun( result$LiqPrem_Tsy_predicted, result$LiqPrem_Tsy )
R2_bond = R2_fun(  result$LiqPrem_bonds_predicted , result$LiqPrem_bonds )
return( data.frame(R2_Tsy=R2_Tsy, R2_bond=R2_bond) )
}
R2s_fun_nested <- function( regs_list,  DT_list,  total_Tsy_IV_vec,  liq_prem_IV_vec ){
R2_Tsy_vec = rep(0, length(regs_list))
R2_bond_vec = rep(0, length(regs_list))
for( iter in 1:length(regs_list) ){
coef_estimations = regs_list[[iter]]$coefficients
result = g_nested(coef_estimations, DT_list[[iter]], prediction=TRUE, total_Tsy_IV=total_Tsy_IV_vec[iter])
# MyPlot( 1:nrow(DT_list[[iter]]), result[,list(LiqPrem_Tsy,LiqPrem_Tsy_predicted)] )
# MyPlot( 1:nrow(DT_list[[iter]]), result[,list(winsorize_down_fun(result$LiqPrem_Tsy_predicted), winsorize_down_fun(result$LiqPrem_Tsy))] )
R2_Tsy_vec[iter] =  summary( lm(winsorize_down_fun(result$LiqPrem_Tsy)~ winsorize_down_fun(result$LiqPrem_Tsy_predicted)) )$r.squared
R2_bond_vec[iter] =  summary( lm(winsorize_down_fun(result$LiqPrem_bonds)~ winsorize_down_fun(result$LiqPrem_bonds_predicted)) )$r.squared
}
return( data.frame( R2_Tsy=R2_Tsy_vec, R2_bond=R2_bond_vec )   )
}
TB_for_reg = TB_preGMM
# note: use the projection of P2CP3M-Deposit rate on FFR as the deposit spread measure.
TB_for_reg[ , CP_deposit_spread:= CP_deposit_spread_projection ]
# Basic plots
par( mfrow=c(1,1) )
DT = TB_for_reg[Year>=1970 , list( P2_CP_Tbill,  CP_deposit_spread , DebtToGDP, Deposits_GDP, VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_P1_spread, DebtToGDP_raw, FFR, YearMon) ]
DT = remove_missing(as.data.frame(DT))
DT1 = TB_for_reg[Year>=1970&Is_Quarter_End , list( P2_CP_Tbill, CP_deposit_spread , DebtToGDP, Deposits_GDP, VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_CP_MMF, DebtToGDP_raw, FFR, YearMon) ]
DT1 = remove_missing(as.data.frame(DT1))
columns = c(1:9)
# N = length(beta_init_list)
counter=0
beta_init_list = list()
for(iter in 1:length(result_core_list)){
for(iter2 in 1:length(result_inner_list)){
if( sum(result_core_list[[iter]]$coefficients>1) + sum(result_inner_list[[iter2]]$coefficients>1) == 0 ){
counter = counter+1
beta_init_list[[counter]] = c( as.numeric(result_core_list[[iter]]$coefficients), as.numeric(result_inner_list[[iter2]]$coefficients) )
}
}
}
regs=list()
regs[[1]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = FALSE) , DT[,columns], beta_init_list,  optfct = "nlminb", type="iterative" )
regs[[2]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = FALSE) , DT1[,columns], beta_init_list, optfct = "nlminb", type="iterative" )
regs[[3]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = TRUE), DT[,columns], beta_init_list, optfct = "nlminb", type="twoStep" )
regs[[4]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = TRUE), DT1[,columns], beta_init_list, optfct = "nlminb", type="twoStep" )
TB_R2 = R2s_fun_nested( regs,  list(DT[,columns],DT1[,columns],DT[,columns],DT1[,columns]), c(FALSE,FALSE,TRUE,TRUE) )
reg_nested_baseline = regs[[4]]
# ===== RFS Revision Table 9 =====
# Output to Latex
stargazer(regs, covariate.labels = c( "$\\rho$", "$\\beta_{\\lambda}$", "$\\eta$", "$\\beta_{\\mu}$" ),
star.cutoffs = c(0,0,0), dep.var.labels.include = FALSE,
dep.var.caption = "Measurement of non-bank liquidity premium",
add.lines = list( c("Total Treasury IVs?", "No","No","Yes","Yes"),
c( "$R^2$ for Tsy Liq Prem", round(TB_R2$R2_Tsy,digits=2) ),
c( "$R^2$ for Q' Liq Prem", round(TB_R2$R2_bond,digits=2) )
),
column.labels=rep(c("P2CP$-$P1CP", "P2CP$-$MMF"),2)
)
# ******** Regressions in "Substitution Between Money and Near-Money Assets". *******
# Created by Wenhao Li, 2021 CC.
# Updated on March 26, 2022.  Ready for RFS Revision Round 2.
# =================== Load Packages ===================
remove(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
source("data_definition_header.R")
# Limit the data period to the end of 2019.
TB = TB[Year<=2019 & !is.na(Deposits_Total)]
# ********** Plots ***********
# ===================== Figures: Time Series Plots for Main Variables =====================
# This figure shows the variations in the explanatory variables in the paper
# ===================== TABLE 1: Log Regressions =====================
index =  !is.element(TB$Year,  WWII_periods) | is.element(TB$Year,  WWII_periods)
index = index & is.element(1:nrow(TB), seq(3,nrow(TB),by=3))   # The second one picks last month in each quarter
# Note: index now is the whole sample
regs=list()
regs[[1]] = lm(  LiqPrem_Log  ~  log(FFR) + log(VIX)  ,  data=TB[index]  )
regs[[2]] = lm(  LiqPrem_Log  ~  log(VIX) + log(DebtToDeposits) ,  data=TB[index]  )
regs[[3]] = lm(  LiqPrem_Log  ~  log(VIX) + log(Debt_to_KVJ_Deposits) ,  data=TB[index]  )
regs[[4]] = lm(  LiqPrem_Log  ~  log(FFR) + log(VIX) + log(DebtToDeposits) ,  data=TB[index]  )
regs[[5]] = lm(  LiqPrem_Log  ~  log(FFR) + log(VIX) + log(Debt_to_KVJ_Deposits) ,  data=TB[index]  )
stargazer( regs, omit = c("f"), omit.stat = c("F", "adj.rsq", "ser"), omit.table.layout="n",
se=RobustSE(regs), df = FALSE, digits = 2 , column.sep.width="0pt", star.cutoffs = c(0,0,0),
covariate.labels = c( "log(FFR$_t$)", "log(VIX$_t $)", "log($ \\frac{\\text{Net Tsy}_t}{\\text{Deposits}_t} $)", "log($ \\frac{\\text{Net Tsy}_t}{\\text{KVJ Deposits}_t} $)"   ),
dep.var.caption = "\\textit{Dependent variable:} log(Treasury liquidity premium$_t$)", dep.var.labels.include = FALSE )
# ===================== TABLE 2: OLS and IV Regressions in Logs =====================
regs_OLS=list()
regs_OLS[[1]] = lm(  LiqPrem_Log  ~  log(FFR)  + log(DebtToDeposits) + log(VIX) ,  data=TB[index]  )
regs_OLS[[2]] = lm(  LiqPrem_Log  ~  log(FFR)  + log(Debt_to_KVJ_Deposits) + log(VIX) ,  data=TB[index]  )
IV_formula_level = "log(FFR) + log(DebtToGDP_raw) + DebtToGDP_raw  + log(VIX)"
regs = list()
regs_monthly_dummies = list()
IV_formulas = list()
quantity_ratios = c("log(DebtToDeposits)",  "log(Debt_to_KVJ_Deposits)" )
for(i in 1:length(quantity_ratios)){
IV_formulas[[i]] = paste0( "LiqPrem_Log ~ log(FFR) + ", quantity_ratios[[i]], " + log(VIX)", " | ",  IV_formula_level)
regs[[i]] = ivreg( IV_formulas[[i]] , data = TB[index])
}
# Output in text
stargazer( c(regs_OLS, regs), type="text", se=c(RobustSE(regs_OLS), RobustSE(regs)),
df = FALSE, digits = 2 , column.sep.width="0pt", star.cutoffs = c(0,0,0) , omit.stat = c("F", "adj.rsq", "ser"),
dep.var.caption = "\\textit{Dependent variable:} log(liquidity premium$_t$)", dep.var.labels.include = FALSE, no.space = TRUE )
## ************************* Instrumental Variable Regressions and Difference Regressions *************************
## This part implements the same instrumental regressions as Nagel for comparison.
# **** Fed Funds Futures Data
TB_FF_Futures = as.data.table(read_xlsx("Fed_Fund_Futures(Bloomberg).xlsx" ))
TB_FF_Futures[ , FF1:=100-`FF1 Comdty` ]
TB_FF_Futures[ , FF2:=100-`FF2 Comdty` ]
TB_FF_Futures[ , FF3:=100-`FF3 Comdty` ]
TB_FF_Futures[ , FF4:=100-`FF4 Comdty` ]
TB_FF_Futures[ , YearMon:=as.yearmon(Date) ]
TB_FF_Futures_monthly = TB_FF_Futures[ ,  lapply(.SD,  mean_fun  )   ,by=YearMon ]    # Note: it is critical how we aggregate the Fed funds futures data into monthly data.
TB_FF_Futures_monthly[ , Date:=NULL ]
TB_FF_Futures_monthly[ ,YearMon := as.Date(YearMon) ]
TB_new = merge( TB,  TB_FF_Futures_monthly, by="YearMon",  all.x=TRUE )
# *** Define the difference variables
TB_Diff = TB_new[ , lapply(.SD, mean_fun) , by=YearMon ]
TB_Diff[ , LiqPrem_Diff:= 100*GenDiff(LiqPrem)   ]
TB_Diff[ , LiqPrem_LogDiff:= GenDiff(LiqPrem_Log)   ]
TB_Diff[ , LiqPrem_LogWinsorized_Diff:= GenLogDiff(LiqPrem_winsorized)   ]
TB_Diff[ , DebtToGDP_LogDiff:= GenLogDiff(DebtToGDP)   ]
TB_Diff[ , DebtToDeposits_LogDiff:= GenLogDiff(DebtToDeposits)   ]
TB_Diff[ , DebtToDeposits_LogDiff_Lag:= GenLogDiff(DebtToDeposits,-1)   ]
TB_Diff[ , FFR_Diff:= GenDiff(FFR)   ]
TB_Diff[ , FFR_LogDiff:= GenLogDiff(FFR)   ]
TB_Diff[ , VIX_Diff:= GenDiff(VIX)   ]
TB_Diff[ , VIX_LogDiff:= GenLogDiff(VIX)   ]
TB_Diff[ , TbillToGDP_LogDiff:= GenLogDiff(TbillToGDP)   ]
TB_Diff[ , TbillToDeposits_LogDiff:= GenLogDiff(TbillToDeposits)   ]
TB_Diff[ , TbillToGDP_LogDiff_Lag:= GenLogDiff(TbillToGDP,-1)   ]  # values a month ago
TB_Diff[ , TbillToDeposits_LogDiff_Lag:= GenLogDiff(TbillToDeposits,-1)   ]  # values a month ago
TB_Diff[ , TbillTo_KVJDeposits:= TbillToGDP/KVJ_fin_sector_debt_to_GDP   ]
TB_Diff[ , TbillTo_KVJDeposits_LogDiff:= GenLogDiff(TbillTo_KVJDeposits)   ]
TB_Diff[ , DebtTo_KVJDeposits:= Debt_to_KVJ_Deposits   ]
TB_Diff[ , DebtTo_KVJDeposits_LogDiff:= GenLogDiff(DebtTo_KVJDeposits)   ]
TB_Diff[ , DebtToGDP_raw_LogDiff := GenLogDiff(DebtToGDP_raw) ]
TB_Diff[ , YearMon:=as.Date(YearMon)]
TB_Diff[ , DummyIndex:=month(as.POSIXlt(YearMon, format="%d/%m/%Y")) ]
TB_Diff[ , FF_implied_change := FF3 - FF2 ]
TB_Diff[ , FF_implied_change := c( NA, FF_implied_change[1:(nrow(TB_Diff)-1)] ) ]
TB_Diff[ , FF_Futures_IV_Level :=  c( NA, FF2[1:(nrow(TB_Diff)-1)] ) ]
# TB_Diff[ , FFfutures_IV := c(NA, FF_implied_change[1:(nrow(TB_Diff)-1)]) ]
TB_Diff[ , FFfutures_IV := FF_implied_change ]
TB_Diff[ , FF_Futures_IV_Level_lag_3M :=  c( rep(NA,3), FF4[1:(nrow(TB_Diff)-3)] ) ]
# *** Regression settings with IVs (monthly dummies)
name_list = paste0( "M", seq(1,11) )
IV_formula = name_list[1]
for( i in 1:length(name_list) ){
TB_Diff[ , (name_list[i]):= DummyIndex==i]
if( i>=2 ){
IV_formula = paste0(IV_formula, " + ", name_list[i])
}
}
IV_formula_monthly_dummies = IV_formula
IV_formula = paste0( IV_formula, " + FFfutures_IV" )
# ===================== TABLE 3: Difference Regressions: OLS and IV ======
# *** Limit to subsamples for comparison with Nagel(2016) on first stage
TB_repoT = TB_Diff[ YearMon >= as.Date("1991-05-01") & YearMon <= as.Date("2020-11-01")  ]  # Note: make the dates fully accurate
index = TB_repoT$YearMon >= as.Date("1991-06-01")   # Note: replicate Nagel's date range.
regs = list()
regs[[1]] = lm( LiqPrem_Diff ~  FFR_Diff + VIX_Diff , data = TB_repoT)
regs[[2]] = lm( LiqPrem_Diff ~  VIX_Diff +  TbillToDeposits_LogDiff , data = TB_repoT)
regs[[3]] = lm( LiqPrem_Diff ~  FFR_Diff + VIX_Diff +  TbillToDeposits_LogDiff , data = TB_repoT)
regs[[4]] = lm( LiqPrem_Diff ~  FFR_Diff + VIX_Diff +  TbillToDeposits_LogDiff + TbillToDeposits_LogDiff_Lag , data = TB_repoT[index])
IV_formula1 = paste0( "LiqPrem_Diff ~ FFR_Diff | ",  IV_formula)
IV_formula2 = paste0( "LiqPrem_Diff ~ TbillToDeposits_LogDiff | ",  IV_formula)
IV_formula3 = paste0( "LiqPrem_Diff ~ FFR_Diff + TbillToDeposits_LogDiff | ",  IV_formula)
IV_formula4 = paste0( "LiqPrem_Diff ~ FFR_Diff + TbillToDeposits_LogDiff + TbillToDeposits_LogDiff_Lag | ",  IV_formula)
regs[[5]] = ivreg( IV_formula1 , data = TB_repoT[index])
regs[[6]] = ivreg( IV_formula2 , data = TB_repoT[index])
regs[[7]] = ivreg( IV_formula3 , data = TB_repoT[index])
regs[[8]] = ivreg( IV_formula4 , data = TB_repoT[index])
stargazer( regs,  se=RobustSE(regs,lag_input = 12) , add.lines = list(  c("Weak instruments test", rep("",8)),  c( "~~CD statistic", CD_stats_vec_from_ivregs(regs) ),  c( "~~Critical value", rep("",4), 6.56 , 6.56, 6.23, 5.87 )  ),   # c("Instruments", rep("",8)), c("~~Futures-implied change", rep("",4), rep("Y",4)), c("~~Seasonal dummies", rep("",4), rep("Y",4)),
df = FALSE, digits = 2 , column.sep.width="0pt", omit.table.layout=c("nd"),   omit.stat =  c("F", "adj.rsq", "ser"), omit="Constant"  , star.cutoffs = c(0,0,0),
covariate.labels = c( "$ \\Delta \\text{FFR}_t $", "$ \\Delta \\log( \\frac{\\text{T-bill}_t}{\\text{Deposits}_t} )  $", "	$ \\Delta \\log( \\frac{\\text{T-bill}_{t-1}}{\\text{Deposits}_{t-1}} ) $", "$ \\Delta \\text{VIX}_t $ " ),
dep.var.caption = "\\textit{Dependent variable: $\\Delta$Treasury liquidity premium$_t$}", no.space=TRUE,
order = c("FFR_Diff", "TbillToDeposits_LogDiff", "TbillToDeposits_LogDiff_Lag", "VIX_Diff" ) )
par(mfrow=c(1,2))
